
### Functions for dynamic simulation (to be loaded before plotting figures)

OneScen.new<-function(obj, ldv, n, scen, sig, num, shocks){
  # CRAN requirements
  times <- sigma.sqr <- alpha.sqr <- NULL
  
  # Create lower and upper percentile bounds of the confidence interval
  Bottom <- (1 - sig)/2
  Top <- 1 - Bottom
  
  # Create data frame to fill in with simulation summaries
  SimSum <- data.frame()
  ShockVals <- data.frame()
  
  # Change data frame for shock values
  for (i in 1:n) {
    if (is.null(shocks))  {
      scenTemp <- scen
    }
    else if (!is.null(shocks)) {
      if (i %in% shocks[, "times"]) {
        scenTemp <- scen
        for (x in names(shocks)[-1]) {
          shocksTemp <- subset(shocks, times == i)
          scenTemp[, x] <- shocksTemp[1, x]
        }
      }
      else if (!(i %in% shocks)) {
        scenTemp <- scen
      }
    }
    
    # Parameter estimates & Variance/Covariance matrix
    Coef <- obj[["coefficients"]][["fixed"]]
    VC <- vcov(obj)
    
    # Draw covariate estimates from the multivariate normal distribution
    Drawn <- mvrnorm(n = num, mu = Coef, Sigma = VC)
    DrawnDF <- data.frame(Drawn)
    
    ##### The intercept
    #names(DrawnDF) <- names(scenTemp)
    
    # Find predicted value
    qiDF <- data.frame(DrawnDF[, 1]) # keep the intercept
    for (u in names(scenTemp)) {
      qiDF[, u] <- DrawnDF[, u] * scenTemp[, u]
    }
    PV <- rowSums(qiDF)
    
    # Create summary data frame
    time <- i
    ldvMean <- mean(PV)
    
    ldvLower <- quantile(PV, prob = Bottom, names = FALSE)
    ldvUpper <- quantile(PV, prob = Top, names = FALSE)
    ldvLower50 <- quantile(PV, prob = 0.25, names = FALSE)
    ldvUpper50 <- quantile(PV, prob = 0.75, names = FALSE)
    
    # Shock variable values
    if (!is.null(shocks)) {
      ShockNames <- names(shocks)[-1]
      TempShock <- scenTemp[, ShockNames]
      ShockVals <- rbind(ShockVals, TempShock)
    }
    # Combine
    TempOut <- cbind(time, ldvMean, ldvLower, ldvUpper, ldvLower50,
                     ldvUpper50)
    SimSum <- rbind(SimSum, TempOut)
    
    # Change lag variable for the next simulation
    scen[, ldv] <- ldvMean
  }
  # Clean up shocks
  if (!is.null(shocks)) {
    CleanNames <- paste0("shock.", ShockNames)
    names(ShockVals) <- CleanNames
    SimSum <- cbind(ShockVals, SimSum)
  }
  col_idx <- grep("time", names(SimSum))
  SimSum <- SimSum[, c(col_idx, (1:ncol(SimSum))[-col_idx])]
  return(SimSum)
}



dynsim.new<-function (obj, ldv, scen, n = 10, sig = 0.95, num = 1000, shocks = NULL, 
                      ...) 
{
  if ("zelig" %in% class(obj)) {
    stop(paste0("dynsim no longer relies on Zelig.\n", "---- Please use `lm`. ----"), 
         call. = FALSE)
  }
  ModCoefNames <- names(obj[["coefficients"]][["fixed"]])
  if (missing(scen)) {
    stop("\nNo scen provided. Please specify scenario values.", 
         call. = FALSE)
  }
  VarMisError <- "\nAt least one variable name in scen was not found in the estimation model."
  if (is.data.frame(scen)) {
    if (any(!(names(scen) %in% ModCoefNames))) {
      stop(VarMisError, call. = FALSE)
    }
  }
  else if (is.list(scen)) {
    for (a in seq_along(scen)) {
      TempDF <- scen[[a]]
      if (any(!(names(TempDF) %in% ModCoefNames))) {
        stop(VarMisError, call. = FALSE)
      }
    }
  }
  if (!is.null(shocks)) {
    if (!is.data.frame(shocks)) {
      stop("\nShocks must be a data frame.", call. = FALSE)
    }
    if (names(shocks)[1] != "times") {
      stop("\nThe first variable of shocks must be called 'times' and contain the shock times.", 
           call. = FALSE)
    }
  }
  if (n <= 0) {
    stop("\nYou must specify at least 1 iteration with the n argument.", 
         call. = FALSE)
  }
  if (sig <= 0 | sig > 1) {
    stop("\nsig must be greater than 0 and not greater than 1.", 
         call. = FALSE)
  }
  if (is.data.frame(scen)) {
    if (nrow(scen) == 1) {
      SimOut <- OneScen.new(obj = obj, ldv = ldv, n = n, num = num, 
                            scen = scen, sig = sig, shocks = shocks)
    }
    else if (nrow(scen > 1)) {
      results <- list()
      for (u in 1:nrow(scen)) {
        scen_sub <- scen[u, ]
        SimTemp <- OneScen.new(obj = obj, ldv = ldv, n = n, 
                               scen = scen_sub, sig = sig, num = num, shocks = shocks)
        results[[u]] <- cbind(scenNumber = u, SimTemp)
      }
      SimOut <- do.call("rbind", results)
    }
  }
  else if (is.list(scen)) {
    results <- list()
    for (u in seq_along(scen)) {
      SimTemp <- OneScen.new(obj = obj, ldv = ldv, n = n, 
                             scen = scen[[u]], sig = sig, num = num, shocks = shocks)
      results[[u]] <- cbind(scenNumber = u, SimTemp)
    }
    SimOut <- do.call("rbind", results)
  }
  class(SimOut) <- c("data.frame", "dynsim")
  return(SimOut)
}
